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Abstract 

This  thesis  is  a  study  of  a  digital  control  technique  known  as 
Output  Predictive  Control  (OPC)  or  Model  Algorithmic  Control  (MAC) .  In 
OPC,  the  behavior  of  the  system  is  predicted  using  its  impulse  response 
function  and  the  desired  response  is  characterized  by  a  reference  tra¬ 
jectory.  Controls  are  computed  iteratively  to  drive  the  system  output 
along  the  desired  trajectory. 

In  an  earlier  study,  the  system  was  made  to  follow  the  reference 
trajectory  exactly,  but  only  at  the  control  application  time;  there  were 
large  oscillations  of  the  output  between  control  changes.  In  this 
study,  the  control  calculation  is  reformulated  as  a  least  squares  curve 
fit  problem,  allowing  some  deviation  from  the  desired  trajectory. 

This  approach  is  applied  as  a  regulator  for  a  very  lightly 
damped  fourth-order  single-input/single-output  system  and  as  a  pitch 
axis  autopilot  in  a  simplified  terrain  following  problem.  A  qualitative 
discussion  of  robustness  properties  is  included. 

The  design  of  the  controller  is  difficult  due  to  the  interrela¬ 
tionships  of  the  internal  parameters;  however,  the  results  of  the 
terrain  following  example  indicate  that  this  is  a  viable  approach  for 
this  problem. 
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THE  APPLICATION  OF  OUTPUT  PREDICTIVE  DIGITAL  CONTROL  TO  WING 


FLUTTER  SUPPRESSION  AND  TERRAIN  FOLLOWING  PROBLEMS 


I  Introduction 


Background 

A  new  control  technique  known  as  Output  Predictive  Control  (OPC) 
or  Model  Algorithmic  Control  (MAC)  has  recently  been  developed  (Ref  6,7, 
8,11,12).  The  control  technique  differs  from  previous  state  variable 
methods  in  that  it  employs  a  discrete  impulse  response  model  rather 
than  a  "state  space  model."  Also,  rather  than  employing  "feedback,"  OPC 
uses  an  explicit  prediction  of  the  future  response,  thus  trying  to  find 
the  future  inj^ut  which  best  matches  a  desired  future  response. 

This  technique  is  conceptually  pleasing  in  that  it  approximates 
a  human  "controller"  in  some  tasks.  As  an  example,  consider  a  pilot 
attempting  to  maintain  a  particular  course.  This  course  can  be  referred 
to  as  a  'set  path'  and  the  control  objective  is  to  maintain  the  set  path. 
If  an  error  exists  between  the  aircraft's  actual  course  and  the  set  path, 
the  pilot  formulates  a  control  input  (or  series  of  control  inputs)  to 
return  the  aircraft  to  the  set  path.  Some  of  the  factors  involved  in 
his  control  calculation  include: 

1.  The  zero-input  response  of  his  aircraft,  i.e.,  the  path  the 
aircraft  will  follow  if  no  additional  controls  are  applied. 
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2.  The  urgency  of  the  error  condition,  i.e.,  how  quickly  must 
he  get  back  to  the  set  path?  The  pilot  will  want  the  error  reduced 
rapidly  if  maintaining  his  present  (incorrect)  course  leads  to  a  hazard¬ 
ous  situation  (i.e.,  collision  with  a  mountain  or  another  aircraft). 

3.  Aircraft  limitations,  his  desire  for  a  smooth  maneuver  or 
passenger  comfort,  and  the  urgency  of  the  situation  will  determine  the 
trajectory  to  fly  from  his  present  position  to  a  point  on  the  set  path. 
The  most  direct  approach  is  to  fly  a  heading  perpendicular  to  the  set 
path.  This  would,  however,  result  in  passing  through  the  set  path  and 
an  even  larger  deviation  or  overshoot  on  the  other  side.  A  more  prac¬ 
tical  solution  is  to  choose  a  trajectory  which  will  reduce  the  course 
error  more  gradually  and  allow  the  pilot  to  roll  out  on  the  set  path. 

4.  The  "system  model"  of  the  aircraft  (i.e.,  how  the  aircraft 
will  react  to  a  given  control  input).  Through  experience,  the  pilot 
knows  his  aircraft's  response  to  a  control  change  at  a  particular  alti¬ 
tude,  airspeed  and  aircraft  configuration. 

Based  on  these  factors,  the  pilot  formulates  a  control  scheme  to  fly  his 
aircraft  to  the  set  path.  After  the  control  has  taken  effect,  he 
assesses  the  results  and,  if  necessary,  reformulates  the  problem. 

During  the  landing  phase,  when  course  (and  glide  path)  control  is  crit¬ 
ical,  the  entire  sequence  can  occur  in  less  than  a  second. 

OPC  or  MAC  is  an  attempt  to  automate  this  type  of  process  using 
a  digital  computer.  The  discrete  impulse  response  model  is  used  both 
to  make  an  explicit  prediction  of  future  output  responses  using  know¬ 
ledge  of  past  inputs  and  to  compute  alternative  future  control  strategies. 
The  one  control  strategy  is  then  selected  which  gives  the  best  match  to 
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a  desired  future  output  trajectory.  The  output  prediction  and  control 
computation  are  performed  closed  loop  on  a  discrete  sample  by  discrete 
sample  basis. 

Original  applications  included  industrial  process  control  where 
a  prii-ri  models  were  not  well  understood  (Ref  12).  In  these  situations, 
on  line  identification  of  the  discrete  impulse  response  function  model 
can  be  employed. 

The  control  strategy  has  also  been  successfully  applied  to  a 
number  of  aerospace  problems  (Ref  6,8)  and  theoretical  investigations 
(Ref  7).  In  many  of  these  aerospace  applications,  the  system  is  very 
lightly  damped,  some  nonminimum  phase.  In  such  situations,  it  was 
found  that  a  direct  application  of  the  output  predictive  control  strat¬ 
egy  to  force  the  output  to  follow  the  desired  path  exactly  would  lead 
to  closed  loop  instability  almost  immediately  (Ref  1).  Specifically, 
it  was  found  that  the  system  could  be  made  to  follow  the  desired  tra¬ 
jectory  exactly,  but  only  at  the  discrete  sample  times  corresponding  to 
a  change  in  control.  The  response  could  be  rapidly  oscillating  between 
control  changes,  leading  to  instability  (Ref  1). 

Problem  and  Scope 

The  objective  of  this  thesis  effort  is  to  formulate  and  test  a 
control  strategy  which  will  not  only  follow  the  desired  trajectory  at 
control  change  times,  but  also  keep  output  oscillations  within  "accept¬ 
able"  bounds  between  control  changes. 
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Single- input/single-output  (SISO)  reduced  order  system  models 
will  be  used  to  investigate  the  control  technique  for  use  as  a  regulator 
for  a  very  lightly  damped  system  and  as  a  pitch  controller  for  a 
fighter-type  aircraft. 

This  report  will  not  go  deeply  into  the  background  theory  of 
OPC;  the  interested  reader  is  referred  to  References  6,  7,  8  and  12. 

This  report  will  focus  on  the  controller  used  to  achieve  the  objective 
as  listed  above. 

A  complete  computer  program  called  IDCOM  (Identification  and 
Command)  has  been  designed  by  the  French  company  of  Andersa/Gerbios  to 
implement  their  version  of  the  output  predictive  technique.  IDCOM  is 
not  available  for  proprietary  reasons,  but  the  controller  implemented 
in  this  thesis  is  conceptually  similar. 

Approach  and  Sequence 
of  Presentation 

The  solution  proposed  in  this  report  is  to  formulate  the  pre¬ 
diction/control  computation  so  that  the  future  output  trajectory  is  cal¬ 
culated  at  a  smaller  sample  spacing  than  the  control  switch  times. 

Trying  to  match  the  desired  path  at  this  finer  spacing  then  creates  an 
overdetermined,  weighted  least  squares  problem  and  allows  prediction 
times  into  the  future  as  long  as  are  computationally  feasible.  The 
internal  energy  states  (causing  oscillations  between  control  changes) 
are  thus  controlled  indirectly  by  keeping  the  future  output  trajectory 
within  a  least  squares  "tube"  rather  than  attempting  to  match  the 
desired  path  exactly  at  only  control  switch  times. 
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The  theory  and  equations  used  in  this  scheme  are  presented  in 
Chapter  II.  Internal  parameters  and  the  effect  of  their  variations  on 
the  sampled  output  and  control  energy  expended  are  discussed  in 
Chapter  III.  A  qualitative  discussion  of  "robustness"  in  the  presence 
of  model  mismatch  is  presented  in  Chipter  IV.  Application  of  the  con¬ 
trol  scheme  as  a  pitch  controller  for  a  reduced  order  model  of  a 
fighter  aircraft  is  presented  in  Chapter  V.  Finally,  conclusions  and 
recommendations  for  further  study  are  listed  in  Chapter  VI. 
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II  Theory 

This  chapter  outlines  the  theory  and  equations  underlying  the 
output  predictive  controller  that  is  developed.  This  development  will 
be  for  the  algorithm  applied  as  a  regulator.  Modifications  for  applica¬ 
tion  of  the  algorithm  to  the  aircraft  pitch  controller  and  terrain 
following  problems  are  discussed  in  Chapter  V. 

Figure  1  is  a  functional  block  diagram  of  the  controller  algo¬ 
rithm.  In  this  formulation  of  OPC,  a  desired  trajectory  (y^)  is  cal¬ 
culated  from  the  point  of  the  'present'  system  output  (y)  to  the  set 
point  <y  )  .  A  control  input  is  calculated  to  zero  out  the  differ¬ 

ence  (z)  botween  the  system's  zero-input  response  (y  and  the 
desired  trajectory  (y^)  .  The  following  sections  explain  each  of  the 

functional  blocks  in  detail. 

Discretization  of  System  Model 

This  thesis  treats  the  discrete  time  control  of  the  nth  order 
linear  time  invariant,  single- input  system  represented  as 

x(t)  =  Ax(t)  +  Bu(t)  (1) 

x(°)  =  Xq 

with  single-output  sampled  at  a  constant  sample  time,  T  , 

y(kT)  =  Cx(kT)  (3) 

Assuming  a  piecewise  constant  input,  the  discrete  time  model  of 


the  system  (1)  -  (3)  is  then  represented  by 


Functional  Block  Diagram  of  the  Output  Predictive  Controller 


w 


where 


x(k+l)  =  Fx(k)  +  Gu(k) 

(A) 

lx 

o 

v-/ 

II 

(5) 

y(k)  =  Cx(k) 

(6) 

FHeAT 

(7) 

G  =  (  /eAT dx)B 

(8) 

The  system  (1)  -  (3)  then  has  the  discrete  impulse  response  sequence 
(Ref  10,11)  denoted  as 

(h(l) ,h(2) ,h(3) ,  ...}  =  (CG,CFG,CF2G,  ...}  (9) 

or 

h(i)  =  CF^G  (10) 

Control  Objectives 

The  controller  is  first  implemented  as  a  regulator.  As  pictured 
in  Figure  2,  the  system  output  is  taken  from  an  arbitrary  initial  point 
(state)  along  a  reference  output  trajectory  to  a  desired  final  value  or 
set  point.  As  an  example,  consider  the  pointing  of  a  radar  antenna 
aboard  a  ship  rolling  and  tossing  In  rough  water.  If  the  desired  antenna 
angle  is  constant  at  @x  (referenced  to  true  north  or  some  inertial 
reference)  the  set  point  for  the  regulator  is  0^  .  The  controller's 

task  then  is  to  drive  the  antenna  to,  and  keep  it  at,  @x  ,  regardless 
of  the  gyrations  of  the  ship. 


5, 

*■ 

.  > 

r 


i 


8 


Fig  2.  Controller  Implemented  as  a  Regulator 


In  the  tracking  application  discussed  in  Chapter  V,  the  control¬ 
ler  is  used  as  a  pitch  axis  autopilot.  Instead  of  a  constant  set  point, 
the  system  is  driven  to  a  time  varying  commanded  pitch  angle  as  dictated 
by  the  requirement  to  fly  along  a  given  terrain  profile. 

Calculate  Points  Along 
Desired  Trajectory 

In  the  regulator  implementation  of  the  proposed  controller,  the 
reference  or  desired  trajectory  is  chosen  as  a  decreasing  exponential 
with  time  constant  Tau.  Discrete  points  along  the  reference  trajectory 
shown  in  Figure  2  can  be  calculated  as 


(11) 


yd(i)  =  yset  '  ^(yset-y) 

where 

y^i)  =  discrete  point  along  desired  path 

yget  =  set  point  or  where  wo  want  the  system  to  go  (for  the 
regulator,  yset=0  ) 

y  =  'present'  output  of  the  system 

=  exp  (  -T^/Tau  ) 

Tj  =  sample  time  at  which  y^(i)  is  calculated 

Tau  =  time  constant  of  first  order  decreasing  exponential 
representing  reference  trajectory  from  y  to  yset 

Figure  2  illustrates  the  situation  for  initial  start-up  of  the  system. 
For  each  subsequent  iteration  of  the  digital  algorithm  the  "present" 
system  output  is  used  for  the  initial  state  and  a  new  reference  trajec¬ 
tory  is  calculated  for  each  iteration. 

Zero-Input  Response 

The  total  system  response  can  be  represented  as 

y(t)  =  yzs(t)  +  yzi(t)  (12) 

where 

yzg(t)  =  zero-state  response 
yzi(t)  =  zero-input  response 

If  the  system  is  started  at  rest  (i.e.,  zero  internal  energy),  the 
response  to  an  input  would  be  termed  the  system's  zero-state  response. 

The  system  response  at  time  k  based  solely  on  past  inputs 
(i.e.,  all  future  inputs  equal  zero),  is  the  zero-input  response.  The 
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zero-input  response  nay  be  found  in  terms  of  the  impulse  response 
function  (Ref  10,11)  or  in  terms  of  the  discrete  system  model  equations 
as 

yzl(i)  =  CF'Sc(O)  (13) 

The  algorithm  is  an  iterative  technique  where  the  "present"  time 
is  always  considered  to  be  t=0  .  When  the  program  is  first  started 

x(0)  is  the  arbitrary  initial  state  of  the  system.  For  each  successive 
iteration  x(0)  is  the  actual  system  state,  thereby  closing  the  loop. 

It  is  beyond  the  scope  of  this  thesis  to  consider  the  "state  estimation 
problem,"  but  clearly  a  closed  loop  estimator  of  the  "state"  is  required 
for  this  particular  implementation  of  the  prediction  calculation. 

Control  Calculation 

It  is  desired  that  the  system  output  match  some  reference  path, 
or  in  equation  form 

yd(i)  =  y(i)  ;  i  =  1,2, ..,n  (14) 

where 

y^(i)  =  points  along  desired  path 
y(i)  =  system  output 

The  discrete  form  of  Eq  (12)  is  substituted  for  the  right  side  of  Eq  (14) 

yd(i)  =  yzi(i)  +  yzs(1)  (15) 

An  expression  relating  yzi(i)  »  the  zero-input  part  of  the  response, 
to  the  current  system  "state"  was  given  by  Eq  (13). 
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Using  the  ideas  of  a  piecewise  constant  input  and  a  discrete 


representation  of  the  convolution  integral  (Ref  10)  gives  an  expression 
for  the  zero-state  response: 


y2S(i) 


i 

l  h(j)u(i-j) 

j=l 


(16) 


Substituting  Eq  (16)  into  Eq  (15)  and  expanding  yields 


yd(l) 

N 

H* 

H-* 

_  ^ _ 1 

h (l)u(0) 

yd(2) 

yzi<2> 

h(2)u(0)+h(l)u(l) 

yd<3) 

• 

= 

y2l(3) 

• 

+ 

h(3)u(0)+h(2)u(l)+h(l)u(2) 

• 

yd(N) 

N 

T.  h(j)u(N-j) 

j-1 

—  m 

(17) 


This  equation  is  solved,  with  no  restrictions,  in  Ref  1  for  the  input  as 


u(0)  = 


yd(1)  "  yzi(1) 
h(  1) 


(18) 


This  is,  in  effect,  one  step  ahead  prediction  and  led  to  instability. 

If  the  problem  is  restricted  to  bring  in  the  idea  of  "smoothing 
terms,"  or  equivalently,  if  the  input  is  held  constant  over  NSM  (number 
of  smoothing  terms)  of  the  desired  points  (  y^(i)  )>  an  overdetermined 

least  squares  problem  evolves.  Figure  3  illustrates  the  concept  of  four 
smoothing  terms  within  the  span  of  one  control  input. 
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A  new  set  of  controls  can  be  defined  as 

[  u(0) ,u(l) , . . ,  u(NSM)  ]  =  u(0) 

[u(NSM+l),u(NSM+-2),..,  u(2*NSM)]  =  u(l) 

•  • 

•  • 

[  .  .  •  ,  u(N)  ]  =  u(L)  (19) 

where  u(0)  is  not  a  row  vector,  but  the  value  assigned  to 
u(0) ,u(l) , . . .u(NSM) . 

Eq  (19)  can  be  substituted  into  Eq  (17),  then  manipulated  to  yield 

z  =  Hu  (20) 
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(23) 

Points  along  the  desired  trajectory,  y^(i)  >  and  the  zero-input 

response  are  still  given  by  Eqs  (11)  and  (13)  respectively,  u  is  a 
vector  of  inputs  with  dimension  "L,"  the  number  of  future  inputs  to 
be  calculated.  H  is  a  matrix  of  impulse  response  functions,  similar  in 
information  content  to  the  Hankel  matrix  (Ref  11). 

Eq  (20)  is  a  linear  equation  which  can  be  solved  for  ju  .  If 
the  control  problem  is  reformulated  using  the  smoothing  terms  approach, 
the  result  is  an  overdetermined  problem;  the  "best"  solution  for  u  is 
a  weighted  least  squares  "curve  fit."  The  cost  function  can  be  defined 
as 


u(0) 

u(l) 


I  u(L) 


J  =  (z-Hu)  'Q (z.~Hu)  (24) 

where  Q  is  a  positive,  semi-definite  weighting  matrix. 

For  computational  simplicity,  one  can  also  choose  for  it  to  be  diagonal. 
The  influence  of  Q  on  the  controller  performance  will  be  discussed  in 
detail  in  Chapter  III. 

One  way  to  solve  the  least  squares  approximation  is  to  develop 
the  Normal  Equation  (Ref  5:122-129;  Ref  9:72,222-224).  Taking  the 
partial  derivative  of  J  with  respect  to  u  ,  setting  the  result  equal 
to  zero  and  solving  for  u  will  yield  the  Normal  Equation: 

.M  ,  -2  (z-iiu)  "QH  =  0  (25) 
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Taking  the  indicated  transpose  and  solving  the  equation  for  u  yields 


z'Qll  =  u'H'QH 

(26) 

u"  =  £>QH(H"QH)“1 

(27) 

u  =  (H'QH)-1H'qz  (28) 


The  dimensions  of  the  elements  of  Eq  (28)  are 
_u  =  An  L  dimension  vector 
H  =  An  (NSM*L)XL  matrix 
Q  =  An  (NSM*L)X(NSM*L)  matrix 
z,  =  An  (NSM*L)  dimension  vector 

Due  to  the  fact  that  Q  is  a  diagonal  matrix  and  all  of  the 
real  information  contained  in  H  is  found  in  the  first  column,  a 
memory  saving  FORTRAN'  coding  technique  was  proposed  by  Dr.  J.  G.  Reid 
and  implemented  by  the  author.  This  technique  is  illustrated  by  the 
solution  of  the  normal  equation  for  a  simple  example  problem  in 
Appendix  A. 

Eq  (28)  is  used  to  solve  for  a  vector  of  inputs  u  .  In  terms 
of  on-line  implementation,  only  the  first  element  of  u  is  used.  The 
problem  is  reformulated  with  each  iteration  and  a  new  u  vector  found. 

Application  of  Control 

After  the  control  has  been  found,  the  system  state  is  updated: 


x(k+l)  =  Fx(k)  +  Gu(0) 


(29 


u(0)  is  the  first  element  of  the  vector  of  controls  found  using  Eq  (28) 
The  output  of  the  system  can  now  be  found  as 


i 

i  Summa:  y 


y(k)  =  Cx(k) 


(30) 


The  general  sample  times  and  index  of  Figure  1  can  now  be 
defined.  The  overall  sample  time,  T  ,  is  the  control  change  time, 

TC  .  The  system  model  is  discretized  at  a  sample  time  T^=TC/NSM 

The  desired  trajectory  (  y^  )  and  zero  input  response  function  (  y  ^  ) 

are  predicted  i=l,2, . . .NSM*L  discrete  points  into  the  future. 

The  relationships  developed  in  this  chapter  have  been  coded  in 
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FORTRAN.  Appendix  B  contains  a  sample  program.  Embodied  in  the  pro¬ 
gram  listed  are  the  options  of  adding  zero-mean  white  Gaussian  noise  to 
the  sampled  input  and/or  control  calculated.  Model  mismatch  options 
are  also  provided.  These  features  are  discussed  in  Chapter  IV. 

A  very  lightly  damped  fourth  order  system  (see  Chapter  III)  is 
chosen  as  a  test  model  for  the  regulator  and  the  results  of  variations 
of  the  internal  parameters  used  by  the  controller 

TC  =  Time  span  of  one  control  input 

Tau  =  Time  constant  of  reference  trajectory 

L  =  Number  of  future  inputs  calculated  per  iteration 

q  —  Weighting  matrix 

NSM  =  Number  of  smoothing  terms  calculated  per  control 

NSM*L  =  Total  number  of  output  points  predicted  into  the 
future 

are  recorded.  The  effect  on  the  system  output  and  control  energy 
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required  to  drive  the  system  to  a  zero  set  point  for  various  combina¬ 
tions  of  internal  parameters  are  the  subjects  of  Chapter  III. 
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Ill  Internal  Parameter  Variation 
and  Results 


The  hypothetical  regulator  described  in  Chapter  II  is  applied 
to  a  lightly  damped  fourth  order  system.  Throughout  this  chapter,  the 
control  objective  is  to  take  the  system  output  from  an  arbitrary 
initial  state,  along  an  exponential  path,  to  a  zero  set  point.  The 
results  of  variations  of  the  internal  parameters  Tau,  Q  ,  NSM  , 

L  and  TC  were  graphically  recorded.  Only  a  few  of  the  infinite  num¬ 
ber  of  possible  combinations  are  presented.  Because  of  the  interde¬ 
pendence  of  the  parameters,  only  a  general  idea  of  the  effects  of 
variations  of  a  particular  parameter  is  possible. 

An  arbitrary  initial  state  and  C  matrix  are  chosen  as 

x^  =  [  20  10  5  0  ]  (31) 

C  =  [  1  0  0  0  )  (32) 

for  demonstration  of  the  characteristics  of  the  controller.  A  random 
number  generator  was  used  at  one  time  to  generate  an  initial  state; 
however,  the  general  characteristics  exhibited  by  the  controller  were 
unchanged . 

System  Model 

An  all  pole  fourth  order  system,  described  by  the  differential 

equation 
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y  (t)  +  1.6  y  (t)  +  274y(t)  +  279y(t)  +  8612y(t)  =  u(t) 


(33) 


is  chosen  as  the  system  model.  y(t)  is  defined  as  an  arbitrary  output 
and  u(t)  as  an  arbitrary  input.  The  eigenvalues 


Xj  2  =  =  “0*25±jl5,4 

(34) 

X^  ^  =  ciju)^  =  -0.55±j6.0 

(35) 

are  representative  of  the  first  and  fifth  aeroelastic  modes  of  the  B-52 
Control  Configured  Vehicle  airspeed  root  locus  (Ref  1:8,18;  Ref  2:2-20). 
Figure  4  shows  the  fourth  order  system's  open  loop  time  response  for  a 
unit  step  input. 

Variation  of  Time  Constant  of  the 
i  Exponential  Reference  Trajectory  - 

Tau 

In  the  regulator  application,  an  exponential  path  is  chosen  as 
the  reference  trajectory  or  desired  path  to  the  set  point.  Tau,  the  time 
constant  of  the  exponential  path,  is  a  measure  of  how  quickly  the  system 
output  is  driven  to  the  set  point. 

Figure  5  shows  the  system  output  for  Tau  =0.1  and  Tau  =  0.3 
Predictably,  as  the  system  output  is  brought  to  the  set  point  more 
rapidly,  the  magnitude  of  the  controls  required  increases.  This  is 
shown  in  Figures  6  and  7. 

Variations  in  the  Weighting 
Matrix  —  Q 

The  Q  matrix  in  Eq  (28)  is  defined  as  a  positive  semi-definite, 
diagonal  weighting  matrix.  In  the  over determined  least  squares  problem 
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Fig  7.  Controls  Applied  for  Tau=0.3 
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OUTPUT  OUTPUT  OUTPUT 


a.  Uniform  Weighting  About  Desired  Trajectory 


TIME 

b.  Closer  Fit  at  End  of  Interval 


TIME 


c.  Closer  Fit  at  Beginning  of  Interval 


Lg  8.  Different  Possible  Output  Constraints  Defined  by 
the  Weighting  Matrix  Q  (Ref  1:79,81-82) 
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formulated  to  find  the  optimal  control  input,  the  relative  sizes  of 
the  elements  of  the  Q  matrix  define  the  shape  of  a  tube  (MIMO)  or 
envelope  (SISO)  along  the  desired  trajectory.  In  the  least  squares 
problem  formulation,  the  output  of  the  system  is  optimized  to  stay  within 
the  envelope  (tube),  as  indicated  in  Figure  8. 

The  weighting  matrix  used  to  generate  the  envelope  illustrated 
by  Figure  8a  is  the  identity  matrix: 


Q3  ~  Diagonal [1*2,4,8,16,...] (nsm*L)X(NSM*L)  (38) 

Q4  =  Diagonal [1,1,.., 1;106,106,..,106](nsmal^x(nsmal)  (39) 

The  progressive  application  of  Q2  ,  Q3  and  Q4  results  in  the  tube 
being  squeezed  tighter  along  the  latter  points  of  the  interval  (i.e.,  the 
oscillation  of  the  system  output  about  the  desired  trajectory  is  being 
progressively  constrained).  For  the  regulator  applied  to  the  fourth 
order  system  discussed  earlier,  the  use  of  Q4  results  in  the  output 
closely  approximating  the  desired  path. 

The  weighting  matrix  denoted  as  Q4  is  used  as  a  "bench  mark" 
for  comparison  of  the  system  output  and  controls  required  for  variations 
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of  the  Q  matrix.  Figures  9  through  17  show  the  system  output  and  con¬ 
trol  inputs  required  for  Q1  ,  Q2  and  Q3  versus  Q4's  use  as  the 
weighting  matrix.  In  the  parameter  identification  boxes  (upper  right 
corner  of  each  graph),  Q1  is  referred  to  as  "Identity  Weighting," 

Q2  is  referred  to  as  "Block  Weighting,"  Q3  as  "Geometric  Weighting" 
and  Q4  as  "MOD  2  Weighting." 

A  desire  to  have  the  system  output  initially  follow  the  desired 
trajectory  very  closely  and  allow  deviations  from  the  desired  path  to 
increase  with  time  is  indicated  in  Figure  8c.  The  weighting  matrix 
formulated  to  do  this  is 

Q5  =  diagonal [...,32,16,8,4,2,1] (nsm*L)X(NSM*L)  (40) 

Figures  18  and  19  show  the  system  output  and  controls  with  Eq  (40)  as 
the  weighting  matrix.  Scaling  of  the  plots  masks  the  fact  that  the  sys¬ 
tem  goes  unstable  almost  immediately. 

The  choice  ef  the  weighting  matrix  to  be  utilized  is,  of  course, 
application  dependent.  Natural  damping  of  the  system,  allowable  peak 
overshoot  and  specified  settling  time  are  all  factors  in  the  selection 
of  Q  .  All  of  these  specifications  are  reflected  in  the  envelope 
(tube)  shape  selected  from  Figure  8.  In  this  regulator  application,  with 
a  very  lightly  damped  system  model,  a  trade-off  between  "nice"  system 
response  and  the  control  energy  required  is  apparent. 

For  the  regulator  applied  to  a  very  lightly  damped  system,  utili¬ 
zation  of  a  weighting  matrix  similar  to  Q3  ,  Eq  (38),  or  Q4  ,  Eq  (39), 
gives  the  best  results.  Use  of  Eq  (40)  is  clearly  a  "bad"  approach. 
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Fig  10.  Controls  Applied  When  Using  Q4  as  the  Weighting  Matrix 


Fig  11.  Controls  Applied  When  Using  Ql  as  the  Weighting  Matrix 
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Fig  12.  Sampled  Output  for  Q4  (A)  and  Q2  Used  as  the  Weighting  Matrix 
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CONTROL  INPUT 


Fig  13.  Controls  Applied  When  Using  Q4  as  the  Weighting  Matrix 
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Controls  Applied  When  Using  Q2  as  the  Weighting  Matrix 


30 


SAMPLED  OUTPUT 
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15.  Sampled  Output  for  Q4  (A)  and  Q3  Used  as  the  Weighting  Matrix 
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Variation  of  Number  of 
Smoothing  Terms  —  NSM 

As  the  number  of  smoothing  terms  per  control  change  is  increased, 
the  system  response  time  increases,  as  seen  in  Figure  20.  The  initial 
increase  from  NSM=2  to  NSM=4  smooths  the  response.  However,  a  fur¬ 
ther  increase  to  NSM=10  reintroduces  oscillations  in  the  output. 

Figures  21  through  23  show  the  controls  applied  for  each  case. 

As  NSM  is  increased,  there  is  a  substantial  decrease  (note  factor  of 
l 

10  on  controls  for  NSM=2  )  in  the  control  energy  required  to  bring 
the  system  to  the  zero  set  point.  As  the  number  of  smoothing  terms 
increases  (i.e. ,  the  reference  trajectory  becomes  better  defined),  less 
control  energy  is  expended  in  controlling  overshoots. 

Variation  of  Number  of  Future 
Controls  Calculated  —  L 

The  solution  of  the  Normal  Equation  (28)  is  a  vector  of  control 
inputs  with  length  or  dimension  L  .  The  quantity  "L*TC"  is  a 
measure  of  how  far  in  the  future  both  prediction  and  calculations  are 
carried  out.  The  vector  of  controls  could  be  sequentially  applied; 
however,  in  this  controller  only  the  first  element  of  the  control  vector 
is  applied  and  the  problem  reformulated  at  each  iteration.  Figure  24 
shows  how  variations  in  L  affect  the  system  output. 

Increasing  from  L=3  to  L=4  damps  the  large  oscillations  and 

h 

increases  the  system  response  tUmq.  Also  investigated  (but  not  shown) 
is  an  increase  to  L=10  which  further  slows  the  system  response  and 
reintroduces  some  low  amplitude  oscillation. 
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Fig  20.  System  Output  for  NSM=*2  (Clean),  NSM=4  (A)  and  NSM=10  (0) 
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SAMPLED  OUTPUT 

rS.OO  -1.00  3.00  7.00  11.00  15.00  19.00  23.00 


Figures  25  through  27  show  that  as  L  is  increased,  there  is  a 
decrease  in  the  control  energy  required  to  bring  the  system  to  the  zero 
set  point. 

The  introduction  of  oscillations  into  the  output  for  L>7  and 
large  reduction  in  control  energy  required  as  L  increases  can  be 
explained  through  an  analysis  of  the  weighting  matrix  (Eq  (38))  uti¬ 
lized,  and  the  fact  that  the  time  interval  in  question  is  increasing  as 
L  increases.  Use  of  Q3  reflects  a  desire  for  a  closer  fit  of  the 
output  to  the  desired  trajectory  toward  the  end  of  the  time  interval. 
This  is  indicated  by  Figure  8b.  As  the  tjime  interval  of  calculation  is 
increased  by  increasing  L  ,  larger  deviations  are  tolerated  early  in 
the  time  interval.  As  less  weighting  is  given  the  early  deviations, 
less  control  energy  is  spent  trying  to  match  the  desired  path  during 
the  period  where  it  (the  desired  path)  is  rapidly  changing. 

A  good  initial  guess  for  L  is  2*N  ,  where  N  is  the  dimen¬ 

sion  of  the  state  vector.  The  final  value  of  L  chosen  depends  on  the 
specifications  for  system  response  time  and  control  input  limitations. 

Variations  in  the  Tir.e  Span  of 
Each  Control  —  TC 

Other  than  a  slight  increase  in  system  response  time,  Figure  28 
indicates  no  real  change  in  system  output  for  variations  in  TC.  There 
is,  however,  a  dramatic  change  in  control  energy  requirements  as  indi¬ 
cated  by  Figures  29  through  31. 

A  direct  design  procedure  for  selection  of  the  "optimal"  sample 
time  to  maximize  closed  loop  robustness  (robustness  is  discussed  in 
Chapter  IV)  is  developed  in  Ref  11.  This  technique  was  designed  for 
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CONTROL  INPUT  *10 
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Fig  29.  Controls  Applied  for  TC=0.075 
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use  with  the  Output  Predictive  Dead-Beat  Controller  (OPDEC)  when 
investigating  model  mismatch.  The  basic  procedure  used  is  to  investi¬ 
gate  the  reciprocal  condition  number,  1/k  ,  of  the  Hankel  matrix 

(H^)  as  a  function  of  sample  time. 


1/k  =  a  .  /a 

min  max 


(41) 


where 

1/k  =  Reciprocal  condition  number 

°min  =  ^n^mum  singular  value  (minimum  magnitude 

T 

eigenvalue  of  H  H  )  of  the  Hankel  matrix 
n  n 

o  «*  Maximum  singular  value  of  the  Hankel  matrix 
max 

Choice  of  the  Hankel  matrix  for  investigation  stems  from 
analyzing  the  equation  used  for  calculation  of  the  inputs  for  the 
OPDEC  (Ref  11): 


V  -  -isi  <«> 

where 

H  =  Hankel  matrix 
n 

u  =  Vector  of  control  inputs 

=  Vector  of  zero-input  responses 

Using  the  same  logic,  the  matrix  in  the  Normal  Equation  to  be  investi¬ 
gated  is  [ii'QH]  . 

The  sample  time  that  maximizes  1/k  is  hypothesized  to  be  the 
"optimal,,  control  change  time  (  TC  )  for  the  minimization  of  pertur¬ 
bation  effects  (model  mismatch,  input  and  measurement  noises)  on  the 
system.  It  is  also  a  realistic  initial  choice  for  TC  ,  without  regard 
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to  any  perturbation  of  the  original  system.  Figure  32  is  a  plot  of 
1/k  of  IH'QH]  for  the  system  model.  Relative  peaks  in  1/k  occur 
at  approximately  0.088  and  0.124  with  valleys  at  0.075  and  0.200. 
Although  Figure  28  shows  no  real  change  in  system  response  for  variance 
of  TC  ,  selection  of  TC=0.124  will  yield  a  much  better  response 
when  the  "original"  system  model  is  in  error  (i.e.,  model  mismatch). 
This  will  be  shown  in  Chapter  IV. 


Nonminimum  Phase  Systems 

A  brief  investigation  of  a  nonminimum  phase  system  (i.e.,  sys¬ 
tem  zeroes  in  the  right-half  plan  (RHP))  is  discussed  in  this  section. 
It  is  hypothesized  that  selection  of  Tau,  Q  and  NSM  is  the  same  as 
for  the  basic  system  with  no  zeroes  in  the  RHP;  however,  because  of  the 
inherent  difficulty  in  controlling  a  nonminimum  phase  system,  the  con¬ 
troller  should  look  further  into  the  future  (increase  in  L  )  for  cal¬ 
culation  of  the  controls.  The  reciprocal  condition  number  of  [H^QH] 
built  from  the  nonminimum  phase  system's  equation  is  used  to  find  TC 

If  a  zero  at  S=+0.3  is  added  to  the  original  system,  the  open 
loop  transfer  function  (OLTF)  becomes 


OLTF  = 


_ (a-0.3) _ 

(s+0.55±j6) (s+0.25±j 15.4) 


(43) 


Figure  33  is  a  plot  of  the  reciprocal  condition  number  of  [hQH]  for 
the  system  model  of  Eq  (43) .  Initially  the  dimension  of  u  is  fixed  at 
L=7  ,  the  "best"  value  for  the  original  (minimum  phase)  system  model. 

Figures  34  through  37  show  the  system  output  and  controls  for  TC=0.125  , 

TC=0.180  (good  sample  times)  and  TC=0.205  (bad  choice,  according  to 
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Fig  32.  Reciprocal  Condition  Number  of  [H^QH]  Using  the 
Basic  System  Model 
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Fig  33.  Reciprocal  Condition  Number  of  [H^QH]  Using  the 
Nonminimum  Phase  System  Model 


47 


CONTROL  INPUT  df  SRMPLPD  OUTPUT 


34.  Sampled  Output  of  Nonminimum  Phase  System  for  TC=0.125 
(Clean),  TC=0.180  (A)  and  TC=0.205  (0);  L=7 


Fig  35.  Controls  for  TC=0.125,  L=7 
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Figure  33).  An  interesting  observation  is  that  for  TC=0.125  ,  a 

supposedly  "good"  sample  time  for  minimization  of  the  effect  of  perturba¬ 
tion  of  the  system  model,  the  system  goes  unstable.  This  is  easier  to 
see  in  Figure  35,  a  plot  of  the  controls  for  TC=0.125  .  Although 

there  is  a  large  transient  with  TC=0.205  ,  the  output  is  driven  to 

(and  kept  at)  the  zero  set  point. 

It  is  hypothesized  that  control  of  the  nonminimum  phase  system 
is  improved  by  allowing  the  controller  to  look  further  into  the  future. 
Figures  38  through  41  show  the  results  of  using  the  same  sample  times 
as  above,  but  increasing  the  dimension  of  iu  to  L=10  .  Instead  of 

instability,  using  TC=0.125  results  in  the  "best"  response  of  the 
three  sample  times. 

This  limited  analysis  has  not  addressed  all  of  the  questions 
concerning  control  of  the  nonminimum  phase  system;  however,  the  main 
factor  in  successful  application  of  the  hypothetical  regulator  to  a 
nonminimum  phase  system  is  an  increase  in  L  over  that  used  to  control 
the  basic  all  pole  system. 

Summary 

An  ironclad  synthesis  technique  is  not  possible  at  this  stage 
of  development  of  the  controller;  however,  the  following  guidelines  are 
presented  for  selection  of  the  internal  parameters  for  control  of  the 
all  pole  system: 

1.  Specifications  for  settling  time,  response  time  and  maximum 
allowable  overshoot  can  be  reflected  in  the  tube  shapes  (Figure  8) 
defined  by  the  elements  of  the  weighting  matrix.  For  the  regulator 
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applied  to  a  lightly  damped  system,  the  weighting  matrices  identified 
in  Eqs  (38)  and  (39)  gives  a  good  system  response.  Eq  (40),  reflecting 
the  tube  shape  indicated  in  Figure  8c  is  clearly  a  bad  choice. 

2.  The  sample  time  which  maximizes  the  reciprocal  condition 
numbe"  of  [H^QH]  is  chosen  as  TC  ,  the  length  of  time  a  control  is 
held  constant.  Although  this  selection  process  may  not  result  in  the 
best  TC  for  an  unperturbed  system  model,  it  is  a  good  initial  choice. 

3.  Tau  is  a  reflection  of  how  quickly  the  error  between  the 
desired  and  actual  paths  is  brought  to  zero.  Selection  of  Tau  is 
dependent  on  the  specified  system  response  time.  Its  value  will  affect 
the  control  energy  expended  in  taking  the  system  output  to  the  desired 
set  point. 

4.  The  number  of  smoothing  terms  (  NSM  )  utilized  depends  on 
the  shape  of  the  desired  path,  A  sufficient  number  of  smoothing  terms 
must  be  utilized  to  characterize  the  reference  trajectory.  For  the 
first  order  exponential  trajectory  utilized  in  the  regulator  application, 
few  smoothing  terms  are  needed.  As  NSM  is  increased,  resulting  in  a 
better  defined  reference  path  for  solution  of  the  least  squares  problem, 
the  controls  needed  to  drive  the  system  output  to  the  set  point  decreases. 
For  the  regulator  example,  2<NSM<7  works  well. 

5.  An  initial  value  for  the  number  of  future  controls  calculated 
(  L  )  per  iteration  is  2n  ;  where  'n'  is  the  order  of  the  system 
model  or  dimension  of  the  state  vector.  The  final  value  of  L  chosen 
has  a  dramatic  effect  on  control  energy  expenditures  as  well  as  system 
response  time.  L=7  is  the  choice  for  the  regulator  example.  Variations 
from  this  introduces  oscillations  in  the  output. 
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For  the  nonminimum  phase  system,  the  same  values  for  Tau,  Q 
and  NSM  are  applicable.  The  sample  time  which  maximizes  the  recipro¬ 
cal  condition  number  of  [H  QH]  for  the  system  equations  is  a  good 
choice  for  TC  .  L  must  be  increased  over  that  used  for  the  all  pole 
system  to  facilitate  control.  Control  of  the  nonminimum  phase  system 
is  not  fully  understood  at  this  point. 

After  an  initial  set  of  internal  variables  is  selected,  a 
scheme  of  parameter  variation  is  necessary  to  get  the  best  set  of  para¬ 
meters  for  the  particular  application. 

There  are  several  interesting  relationships  between  this 
"smoothing  approach"  and  the  Output  Predictive  Dead-beat  Controller 
(OPDEC)  discussed  in  Refs  3  and  11.  Selection  of  L=2n  and  Q=Q4 
given  by  Eq  (39)  for  implementation  of  the  "smoothing  approach"  is 
actually  an  approximation  of  the  scheme  used  in  the  OPDEC  approach  and 
gives  much  the  same  kind  of  performance  as  OPDEC  itself. 
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IV  Perturbations  of  the  System  Model 

For  many  practical  problems,  the  true  system  model  is  unknown 
or  the  system  dynamics  change  with  environmental  conditions.  An 
example  of  such  a  system  is  an  aircraft  for  which  the  flight  character¬ 
istics  change  as  a  function  of  altitude,  airspeed,  attitude  or  some  other 
parameter.  Either  of  the  above  situations  can  be  considered  a  case  of 
model  mismatch;  that  is,  the  actual  system  equations  differ  from  the 
system  model  used  in  designing  the  controller.  In  this  thesis,  model 
mismatch  is  simulated  by  perturbation  of  the  eigenvalues  of  the  plant 
or  "A"  matrix.  The  technique  for  selection  of  the  best  value  of  TC 
(see  Chapter  III)  for  closed  loop  "robustness"  is  tested  in  the  model 
mismatch  simulation. 

The  original  example  system  model  developed  in  Chapter  III  is 
used  in  this  chapter.  The  controller  is  applied  as  a  regulator  with 
the  control  objectives  as  listed  in  Chapter  II.  The  initial  state  and 
output  relationship  are  as  given  by  Eqs  (31)  and  (32). 

Robustness 

A  robust  controller  is  one  that  continues  to  perform  properly 
despite  system  perturbations  (i.e.,  model  mismatch  and  noise  corrupted 
information  fed  back  to  the  controller) .  For  3n  analytical  discussion 
of  robustness,  see  Refs  7  and  11.  For  this  analysis  robustness  is 
measured  in  terms  of  stability  only. 
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Model  Mismatch 


For  discussion  of  model  mismatch,  several  system  models  are 
defined  -  the  Truth  model  and  Delta  models.  The  Truth  model  is  the 
original  system  represented  by  Eq  (33) .  A  Delta  model  is  formed  by 
changing  the  natural  frequencies  associated  with  each  of  the  eigenvalues 
of  the  original  system  by  10%,  20%  or  30%.  Table  I  lists  the  OLTFs  for 
each  model. 


TABLE  I 

System  Models  Utilized  in  Model  Mismatch  Analysis 


System  Model 


OLTF 


Truth 


1 


(s+C.55±j6)  (s+0.25±j 15 .4) 


Delta  1 

(+10%  eigenvalue  change) 


1 


(s+0. 605±j6 .6)  (s+0.275±j 16.94) 


Delta  2 

(+20%  eigenvalue  change) 


1 


(s+0.66±j7.2)  (s+0. 3±j 18.4) 


Delta  3 

(+30%  eigenvalue  change) 


1 


(s+0. 715±j 7.8)  (s+0. 325±j20. 2) 


The  regulator  is,  of  course,  designed  for  use  with  the  Truth 
model.  Model  mismatch  occurs  whenever  the  discrete  matrix  representa¬ 
tion  of  one  of  the  Delta  models  is  used  in  calculation  of  the  control 
and/or  prediction  of  the  zero-input  response.  One  entry  point  into  the 
control  calculation  is  through  the  H  matrix,  given  by  Eq  (22).  Each 
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element  of  H  is  given  by 


h(i)  =  CF±-1G  (10) 

The  zero-input  response  prediction  is  given  by 

yzl(i)  =  CF^x(O)  (13) 

Both  H  and  y  ^  are  contained  within  the  Normal  Equation  (28)  used  to 
find  the  input.  In  this  analysis,  however,  they  are  considered  sepa¬ 
rately  as  independent  entry  points  for  the  discrete  matrices  describing 
the  Delta  models. 

Model  Mismatch  Simulation 
Results 

The  basic  regulator  described  in  Chapter  II  is  augmented  with  a 
series  of  logical  switches  (see  Appendix  A)  allowing  the  Delta  model 
matrices  to  be  used  for  prediction  of  yzi  and/or  calculation  of  H  . 
Tables  II  and  III  summarize  the  results  of  the  simulation  for  selected 
values  of  TC  .  Figures  showing  the  system  output  and  controls  applied 
are  in  Appendix  C.  Figures  are  included  for  those  cases  where  the 
system  is  slowly  going  unstable. 

Selection  of  the  various  values  of  TC  for  the  model  mismatch 
simulation  was  made  with  Figure  42  in  mind.  In  Chapter  III  it  is 
hypothesized  that  the  discrete  time  which  maximizes  1/k  of  [h  QH]  is 
the  best  TC  for  closed  loop  robustness  of  the  controller.  In  relation 
to  the  value  of  TC  ,  the  following  observations  are  made  concerning 


Tables  II  and  III: 


TABLE  II 


Summary  of  Results  for  Perturbed  or  Delta  Models  Used  to  Calculate 

the  Zero-Input  Response 


System  Model 

TC  (sec) 

Results 

Delta  1 

| 

j 

0.075 

stable  (Figure  C-l) 

(10%  eigenvalue  change) 

0.088 

stable  (Figure  C-4) 

0.124 

stable  (Figure  C-7) 

0.200 

unstable 

Delta  2 

0.075 

stable  (Figure  C-10) 

(20%  eigenvalue  change) 

0.088 

stable  (Figure  C-l 3) 

0.124 

unstable  (Figure  C-16) 

0.200 

unstable 

Delta  3 

0.075 

stable  (Figure  C-19) 

(30%  eigenvalue  change) 

0.088 

stable  (Figure  C-22) 

0.124 

unstable 

unstable 


TABLE  III 


Summary  of  Results  for  Perturbed  or  Delta  Models 
Used  to  Find  H 


System  Model  TC  (sec)  Results 


Delta  1 

(10%  eigenvalue  change) 

0.075 

0.088 

0.124 

0.200 

stable  (Figure  C-25) 

stable  (Figure  C-28) 

stable  (Figure  C-31) 

unstable 

Delta  2 

0.075 

unstable  (Figure  C-34) 

(20%  eigenvalue  change) 

0.088 

stable  (Figure  C-37) 

0.124 

stable  (Figure  C-40) 

0.200 

unstable 

Delta  3 

0.075 

unstable  (Figure  C-43) 

(30%  eigenvalue  change) 

0.088 

stable  (Figure  C-46) 

0.124 

stable  (Figure  C-49) 

0.200 

unstable 

00  '<f 


Fig  42.  Reciprocal  Condition  Number  of  [H^QH]  for  the 
"Truth"  or  Basic  System  Model 
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1.  From  Table  II,  stability  is  maintained  for  all  Delta  models 
only  with  the  smallest  values  of  TC  (0.075  and  0.088). 

2.  For  TC=0.200  ,  the  system  is  unstable  for  all  Delta 

models  to  both  H  and  y 

3.  From  Table  III,  the  only  control  change  times  for  which 
stability  is  maintained  for  all  Delta  models  correspond  to  the  "peaks" 
in  Figure  42. 

As  TC  increases,  the  time  span  of  the  zero  input  response 

prediction  also  increases,  thus  amplifying  the  error  involved  when 

using  a  Delta  model  for  prediction  of  y  .  .  For  any  practical  appli- 

Z  X 

cation  of  the  controller,  closed  loop  prediction  would  be  accomplished 
through  implementation  of  an  observer  or  Kalman  filter  for  use  directly 
in  the  control  law  (17)  or  (28).  This  involves  a  transformation  to  an 
"output  predictive  state  coordinate  system,"  an  idea  developed  by 
Dr.  Reid  in  Ref  11.  The  new  "state  vector”  is  actually  the  predicted 
output  (see  Ref  11).  When  the  zero  input  prediction  responsibility  is 
assumed  by  a  Kalman  filter,  the  results  of  Table  II  will  be  obviated. 
Selection  of  TC  can  be  based  on  the  information  contained  in  Figure  42, 
substantiated  by  the  results  of  Table  III. 

At  TC=0.200,  1/k  is  very  small  (see  Figure  42),  reflecting  the 
fact  that  [H^QH]  is  ill-conditioned  in  terms  of  robustness.  This  is 
confirmed  by  Table  III,  in  that  for  TC=0.200  the  system  is  unstable 
for  all  Delta  models  sent  to  the  calculation  of  H 

TC=0.088  and  TC=0.124  correspond  to  the  two  largest  peaks 
in  the  reciprocal  condition  number,  Figure  42.  As  seen  in  Table  III, 
with  these  values  of  TC  the  system  maintains  stability  with  up  to  a 
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30%  change  in  eigenvalues.  Further,  comparison  of  the  system  response 
for  these  two  values  of  TC  (Figures  43  and  44)  indicates  that  the 
controller  is  "more  robust"  at  the  time  which  maximizes  1/x  (TC=0.124) 


Summary 

The  model  mismatch  simulation  points  out  two  important  factors 
to  be  considered  in  a  "real  world"  application.  The  first  is  that  to 
avoid  the  results  indicated  in  Table  II,  a  Kalman  filter  or  observer 
should  be  implemented  to  perform  the  future  state  estimation  function. 
Second,  in  applications  calling  for  a  robust  controller,  the  selection 
of  control  change  tine  (  TC  )  is  a  critical  design  step.  The  best 
choice  for  TC  is  the  discrete  time  which  maximizes  the  reciprocal 
condition  number  of  [H^QH]  .  The  results  demonstrated  in  this 
chapter  concerning  the  selection  of  TC  corresponds  to  the  same  result 
derived  analytically  in  Ref  11  for  the  case  of  OPDEC. 

This  chapter  concludes  the  discussion  of  the  regulator  applica¬ 
tion  of  the  "smoothing  approach"  to  Output  Predictive  Control. 

Chapter  V  discusses  the  technique  implemented  as  a  pitch  controller  on 
a  modern  aircraft  in  a  terrain  following  problem. 


V  Controller  Algorithm  Implemented  as  a 
Pitch  Axis  Autopilot 

In  this  chapter  the  algorithm  developed  in  Chapter  II  is  modi¬ 
fied  : lightly  for  implementation  as  a  pitch  axis  autopilot  for  a  terrain 
following  problem.  The  control  objectives  and  aircraft  model  are 
discussed  first,  followed  by  an  explanation  of  the  changes  necessary  to 
the  algorithm  developed  in  previous  chapters.  Two  test  cases  are  pre¬ 
sented  and  the  results  discussed.  Model  mismatch  was  not  investigated. 
This  development  is  not  an  attempt  to  present  an  ideal  solution  to  the 
terrain  following  problem,  but  is  included  in  the  thesis  as  an  initial 
investigation  into  a  possible  application  area  for  the  smoothing 
approach  to  Output  Predictive  Control. 

/ 

Control  Objectives 

In  Chapters  III  and  IV  the  smoothing  approach  to  Output  Predic¬ 
tive  Control  is  implemented  as  a  regulator.  In  that  case  the  set  point 
is  constant  at  zero  for  all  time.  In  this  chapter,  the  controller  is 
used  as  a  pitch  axis  autopilot  and  the  objective  is  to  follow  a  desired 
altitude  profile;  the  set  point  becomes  a  time  varying  "set  path." 

System  Model 

The  equations  of  motion  for  an  aircraft  are  nonlinear  with  time 
varying  coefficients.  By  assuming  small  perturbations  about  a  specific 
operating  point,  however,  a  linearized  set  of  equations  is  obtained. 

A  linearized  version  of  the  longitudinal  dynamics  (called  the  short 
period  approximation)  for  a  modern  fighter  aircraft,  linearized  about  a 
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straight  and  level  operating  point  at  0.8  Mach,  10,000  feet  altitude 


and  Uo=constant=1077  fps  ,  is  given  by 


q 
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where 


(44) 


q  =  Pitch  rate  (radians/sec) 

“  =  Angle  of  attack  (radians) 

0  =  Pitch  angle  (radians) 

£  =  Elevator  deflection  (radians) 

The  elevator  deflection  is  assumed  limited  to  -0.262  rads£££0.426  rads 
For  the  short  period  approximation,  f=0.69  and  wn=4.8  rad/sec. 
Figure  45  shows  the  angles  listed  above  and  includes  the  flight  path 
angle  (  y  ),  the  angle  between  the  horizon  and  flight  path. 


Fig  45.  Perturbation  Angles  for  Linearized  Longitudinal  Dynamics 
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With  the  forward  velocity  (  UQ  assumed  constant  in  short 
period  approximation)  maintained  constant  through  use  of  an  automatic 
throttle  system  or  pilot  inputs,  the  additional  state  of  altitude 
(  h  )  is  added  through  the  relationship 


h(t)  =  UQY(t)  =  uo[0(t)-«(t>]  fps 


(45) 


where 


h(t)  =  Vertical  velocity  (fps) 

U  =  Constant  =  1077  fps 
o  v 

The  system  model  becomes 
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(46) 


The  transfer  function  of  interest  is  altitude  change  for  a  given 
elevator  deflection 


h(s)  _  -356 .81 (s-14 ■ 93)  (s+18.78) 

E(S)  sz(s+3. 34±j3.46)  (47) 

The  controller,  then,  will  use  Eq  (46)  as  the  system  model  to  calculate 
the  elevator  deflection  necessary  to  follow  a  desired  altitude  profile. 

Implementation 

The  algorithm  used  for  the  regulator  is  modified  slightly  for 
the  terrain  following  problem.  These  changes  include  the  generation  of 
a  set  path,  a  modification  of  how  the  system  error  is  reduced,  the 
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NOMINAL  ALT  I TUOE ( FT )  «10 


TABLE  IV 


Effective  Size  of  Obstacle  for  Variation  of  the  Frequency 
(ui)  of  the  Sinusoidal  Set  Path 


Frequency  of  Oscillation 
(to)  Hz 

Base  of  Obstacle 
(.t) 

Height 

(Ft) 

0.50 

6772 

1000 

0.60 

5643 

1000 

0.75 

4515 

1000 

0.90 

3762 

1000 

1.00 

3386 

1000 

Fig  47.  Set  Path  or  Desired  Altitude  Profile  for  A=1000  ft, 
w=0.5Hz,  and  TC=0.5  sec 
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Modification  of  System  Error  Reduction  Relationship.  In 


Chapter  II,  Eq  (11)  is  given  as  the  relationship  for  calculation  of 
points  along  the  desired  trajectory  from  present  position  (  y  )  to 
the  set  point  (  yget=0  for  regulator) : 

yd(1)  =  yset-t*i(yset"y)  ;  1  =  1»2>---NSM*L  (ID 

Eq  (11)  is  based  on  the  fact  that,  for  the  regulator  application,  the 
set  point  is  constant  for  all  time.  A  modification  of  this  relationship 
is  necessary  because  there  is  no  longer  a  constant  set  point,  but  a  time 
varying  set  path. 

Assuming  that  an  error  exists  between  the  actual  system  output 
and  the  set  path,  then  y^(t)  is  a  function  of  the  present  system  out¬ 
put,  present  set  point  along  the  set  path,  future  set  points  and  how  we 
wish  to  reduce  the  system  error.  If  the  error  is  to  be  exponentially 
reduced,  then 


Desired  Future  Error  =  3  (Present  Error)  (50) 

where  3  is  as  given  in  Eq  (11).  Substituting  into  Eq  (49) 

y  ,(t  +iT.)-v  .  (t  +iT-.)  =  ty ( t  )~y  .  (t  )]  (51) 

'dv  now  1'  -set  now  1  now  'set  now' 

F>r  the  case  of  a  constant  set  point  where  y  (t  +iT)  =  y  (t  )  , 

be  L  IK/W  be  L  11UW 

Eq  (52)  reduces  to  Eq  (11).  Eq  (51)  is  used  for  calculation  of  discrete 
points  along  the  desired  trajectory. 

Selection  of  Internal  Parameters.  As  discussed  in  Chapter  III, 
tin  of  the  Internal  parameters  is  dependent  on  the  physical 
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system.  In  this  application  of  the  algorithm,  the  parameters  are 
chosen  with  physical  restrictions  in  mind. 

In  an  effort  to  use  a  physically  realizable  control  change  time, 
TC  is  chosen  as  TC=0.5 

Tau,  the  time  constant  of  the  decreasing  exponential  of  Eqs  (11) 
and  (50) ,  is  a  measure  of  how  quickly  the  error  between  the  desired 
and  actual  altitude  profile  is  reduced.  For  this  example,  Tau  is 
chosen  as  Tau=0.3  .  As  Tau  is  decreased  from  this  number,  aircraft 

accelerations  become  excessive. 

A  sufficient  number  of  smoothing  terms  (  NSM  )  to  give  a  well 
defined  characterization  of  the  desired  or  set  path  must  be  utilized. 
Values  over  the  range  3<NSM<8  were  tried,  with  NSM=4  chosen.  As 
the  frequency  of  oscillation  of  the  sine  wave  used  as  the  set  path 
increases,  more  smoothing  terms  per  control  are  required  to  characterize 
this  path. 

Values  of  L  ,  the  number  of  future  control  changes  calculated 
per  iteration,  over  the  range  4<L<10  were  tried.  A  value  of  L=7  is 
chosen  to  effect  a  good  tradeoff  between  response  time  and  the  magnitude 
of  the  elevator  deflections  calculated. 

The  choice  of  a  weighting  matrix  for  the  Normal  Equation  (28) 
is,  of  course,  dependent  on  the  values  of  L  and  TC  .  Each  of  the 
weighting  matrices  discussed  in  Chapter  III  and  Eq  (52)  were  tested  as 
Q  .  Q6  is  basically  Q2  [Eq  (37)]  with  the  column  order  reversed 
and  the  nonzero  elements  moved  to  the  principal  diagonal: 

Q6  =  Diagonal [64,64,64,64 ;32,32, 32,32 jl6, ... ,l]2gX28  (52) 
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where  each  block  Is  NSM=4  units  long.  The  resulting  average  and 
peak  altitude  errors  for  each  weighting  matrix  are  listed  in  Table  ,V. 
Based  on  this  data  Q3  ,  Eq  (38),  is  chosen  as  the  weighting  matrix 
to  be  utilized.  Graphical  results  for  the  case  of  Q=Q4  and  w=0.5 
are  shown  in  Figures  48  through  49.  It  is  interesting  to  note  that 
although  Q4  was  the  best  choice  for  the  regulator  application,  its 
use  induces  large  altitude  errors  in  this  example. 

Final  Tests  and  Results 

A  sum  of  sinusoids  is  chosen  as  the  set  path  or  desired  altitude 
profile  for  final  testing  of  the  controller: 

y  (i)  =  Asinii)j(iTj)  +  Asin^CiT^)  +  Asinw^(iTj)  (53) 

The  results  of  two  test  cases  are  presented.  For  both  tests, 
A=1000  feet  and  the  aircraft  is  started  from  an  initial  condition  of 
straight  and  level  flight.  For  the  first  case  0)^=0. 5  ,  O)2=0.6  and 

0)^=0. 2  ;  for  the  second,  0)^  and  U)2  are  unchanged  and  to^O.S 

Table  VI  summarizes  the  results  of  the  tests. 

Figures  50  through  57  show  the  results  of  the  test  cases. 

Included  for  each  test  are  plots  of  the  desired  and  achieved  altitude 
profiles,  controls  applied  and  pitch  angle  variations.  Plots  of  the 
achieved  altitude  and  control  inputs  show  that  the  controller  came  to 
"steady  state"  tracking  long  before  the  inputs  become  steady  state; 
this  is  considered  a  significant  result. 

Examination  of  the  pitch  angle  history  (Figures  53  and  57)  for 
either  case  reveals  angles  well  outside  the  range  of  "small  perturbation" 
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TABLE  V 


Average  and  Peak  Altitude  Error  for  Variation  of  the  Frequency  of 
Oscillation  (u)  of  the  Sinusoidal  Set  Path  and  Weighting 

Matrix  Utilized 


Weighting  Matrix 

to(Hz) 

Peak  Error  (ft) 

Average  Error  (ft) 

Q1 

mm 

66.42 

35.61 

K  mi 

79.49 

40.65 

Eq  (36) 

131.09 

51.39 

■M 

263.81 

74.42 

1.00 

1051.37 

160.80 

Q2 

mm 

68.17 

35.38 

81.38 

40.43 

Eq  (37) 

164.04 

51.84 

378.67 

83.44 

1.00 

1062.53 

159.80 

Q3 

mm 

66.51 

36.75 

mSm  ■ 

76.96 

42.35 

Eq  (38) 

96.05 

52.55 

mSSM 

113.70 

61.98 

1.00 

124.23 

68.86 

Q4 

mm 

1007.79 

100.04 

■SHR 

1007.79 

104.12 

Eq  (39) 

mBEEm 

1007.79 

114.60 

HI' 

1007.79 

126.70 

1.00 

1007.79 

136.03 

Q5 

mm 

333.72 

74.52 

213.09 

59.75 

Eq  (40) 

Wmm 

1406.44 

214.43 

WKXM 

1116.40 

368.62 

1.00 

6421.20 

1766.51 

Q6 

mm 

67.52 

36.01 

bob 

81.02 

41.10 

Eq  (52) 

mbm 

101.15 

51.09 

162.67 

62.93 

1.00 

703.38 

161.68 

Fig  48.  Set  Altitude  (A)  and  Achieved  Altitude  (Clean)  for 
Q=Q4  and  U)=0.5  Hz  and  TC=0.5  sec 


CONTROL  INCUTS 
TC«.5  TftUx.S  NSttM 
L»7  Ub.S  IXIOOO 
I.C.MTbLTICL  0>04 


’o.oo  4.00  8.00  12.00  16.00  20.00  24.00  28.00 

TIME  IN  SECONDS 


Fig  49.  Control  Inputs  for  Q=Q4  and  (0*0.5  Hz  and  TC=0.5  sec 
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TABLE  VI 


Sinusoidal  Frequencies  Used  and  Peak  and  Average 
Altitude  Errors  for  the  Test  Cases 


Variable 

Test  Case  1 

Test  Case  2 

A(f  t) 

1000 

1000 

(^(Hz) 

0.5 

0.5 

u>2(Hz) 

0.6 

0.6 

o>3(Hz) 

0.2 

0.8 

Average  Error  (ft) 

69.31 

82.18 

Peak  Error  (ft) 

171.59 

254.96 

Fig  50.  Set  Path  or  Desired  Altitude  Profile  for  Test 
Case  1  (see  Table  VI) 


Fig  51.  Desired  Altitude  (A)  and  Achieved  Altitude  (Clean) 
for  Test  Case  1  (see  Table  VI) 


75 


Fig  52.  Control  Inputs  for  Test  Case  1  (see  Table  VI) 


Fig  53.  Pitch  Angle  (9)  for  Test  Case  1  (see  Table  VI) 
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Fig  54.  Set  Path  or  Desired  Altitude  Profile  for  Test  Case  2 
(see  Table  VI) 


Fig  55.  Desired  Altitude  (A)  and  Achieved  Altitude  (Clean) 
for  Test  Case  2  (see  Table  VI) 
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Fig  56.  Control  Inputs  for  Test  Case  2  (see  Table  VI) 


Fig  57.  Pitch  Angle  (0)  for  Test  Case  2  (see  Table  VI) 
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angle  for  which  the  system  model  is  valid.  Although  the  altitude 
errors  indicated  in  Figures  51  or  55  cannot  be  used  to  show  that  the 
controller/aircraft  could  actually  fly  the  test  profiles  as  indicated, 
the  test  results  indicate  that  this  controller  algorithm  is  a  viable 
one  for  the  terrain  following  problem. 

It  is  concluded  that  the  algorithm  is  readily  adaptable  to  a 
tracking  problem.  Limitations  due  to  the  linearized  system  model, 
however,  make  the  results  of  the  terrain  following  tests  unusable 
except  for  comparison  purposes.  In  order  to  get  more  meaningful  results 
the  problem  should  be  reformulated  to  run  altitude  as  a  function  of 
down-range  position  and  utilization  of  one  of  the  following  alternatives 

1.  Go  to  on-line  identification  of  the  impulse  function  and 
use  of  a  Kalman  Filter  for  prediction  purposes. 

2.  A  series  of  linearized  models  can  be  stored,  each  to  be 
used  within  the  appropriate  range  of  flight  conditions  (similar  in  con¬ 
cept  to  gain  scheduling) . 

3.  Make  the  controller  so  "robust"  that  the  model  mismatch 
has  little  effect. 

Appendix  D  contains  the  FORTRAN  code  used  for  Test  Case  1. 


VI  Conclusions  and  Recommendations 


Chapter  VI  contains  the  conclusions  drawn  from  this  study  of 
the  smoothing  approach  to  Output  Predictive  Control  and  some  possible 
areas  for  further  study. 

Conclusions 

In  Ref  1,  the  control  problem  is  formulated  for  exact  matching 
(at  control  change  time)  of  the  system  output  to  the  desired  trajectory 
and  prediction  of  the  2ero  input  response  is  for  one  step  ahead.  For 
a  system  higher  than  second  order  the  above  formulation  led  to  insta¬ 
bility.  This  problem  was  corrected  by  looking  further  into  the  future 
for  control  calculation  and  by  allowing  some  deviation  of  the  output 
around  the  desired  trajectory. 

For  the  smoothing  approach  implemented  as  a  regulator  the  results 
are  good  -  even  in  the  face  of  30%  model  mismatch.  The  level  of  success 
or  goodness,  however,  is  dependent  on  the  proper  selection  of  the 
internal  parameters  of  the  regulator.  The  "best"  set  of  parameters 
for  the  smoothing  approach  implementation  approximates  the  Output  Pre¬ 
dictive  Dead-Beat  Controller  (OPDEC)  found  in  Refs  3  and  11.  The  most 
attractive  feature  of  OPDEC  is  that  there  is  only  one  internal  parameter 
to  be  chosen,  the  control  change  time  (TC) ,  and  a  direct  design  pro¬ 
cedure  has  been  developed  for  the  selection  of  TC 

Investigation  of  the  terrain  following  application  indicates 
that  the  algorithm  is  easily  adapted  to  a  tracking  task.  Limitations 
due  to  the  linearized  system  model,  however,  make  the  current  results 
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of  the  example  inconclusive  for  a  practical  aircraft  controller  with¬ 
out  the  probable  use  of  model  adaption.  However,  this  example  illus¬ 
trates  good  performance  of  the  output  predictive  controller  in  a 
demanding  tracking  situation. 

Recommendations  for  Further 
Study 

Investigation  in  the  following  areas  could  provide  further 
insights  into  the  general  area  of  Output  Predictive  Control: 

1.  Further  development  of  the  synthesis  technique  for  selection 
of  the  internal  parameters  of  the  controller  (i.e.,  TC  ,  L  ,  NSM  , 
Q  and  Tau  ) .  A  clearer  understanding  of  the  interrelationships  of 
the  parameters  will  facilitate  application  of  the  controller  to  new 
problems . 

2.  An  indepth  study  of  the  robustness  properties  of  the 
algorithm  is  necessary  for  quantification  of  these  properties. 

3.  A  reformulation  of  the  tracking  problem  for  more  meaningful 
results  would  include  running  altitude  as  a  function  of  downrange 
position  and  utilization  of  one  of  the  following  alternatives: 

a.  Implementation  of  on-line  identification  of  the  impulse 
response  function  and  use  of  a  Kalman  Filter  for  prediction. 

b.  A  series  of  linearized  models  can  be  stored,  each  to  be 
used  within  the  appropriate  range  of  flight  conditions  (similar  in 
concept  to  gain  scheduling) . 

c.  Insure  that  the  controller  is  so  robust  that  the  model 
mismatch  introduced  by  using  a  linearized  model  outside  its  range 

has  little  effect  on  the  closed  loop  performance. 
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Appendix  A 


Procedure  for  Solution  of  the  Normal  Equation 


The  procedure  for  solving  the  Normal  Equation  for  a  simple 
example  problem  is  presented  to  illustrate  the  dimensions  of  its  compo¬ 
nents  and  to  show  a  series  of  matrix  manipulations  which  result  in  a 
savings  of  computer  storage  space.  Eq  (28),  the  Normal  Equation,  is 
repeated  for  convenience: 


u  =  (H"QH)_1H'Qz 


(28) 


where 


u  =  An  L  dimension  vector 
H  =  An  (  NSM*L  )  X  L  matrix 

Q  =  An  (  NSM*L  )  X  (  NSM*L  )  diagonal  matrix 


*  =  ^d  "  *zi 


(21) 


£  =  An  (  NSM*L  )  dimension  vector 

If  three  future  control  inputs  per  iteration  are  calculated 
(  L=3  ) ,  three  smoothing  terms  per  control  input  are  used  (  NSM=3  ) 
and  Q  is  chosen  as  Q3  ,  Eq  (38) ,  a  weighting  matrix  discussed  in 
Chapter  III,  then 


8A 


"hCD  | 

0 

;  0 

h(l)+h(2)  ' 

0 

1  0 

h(l)+h(2)+h(3)| 

0 

'  0 

•  * 

1 

h(l) 

j  0 

•  l 

b(l)+ h(2) 

«  0 

6  1 

l  h(i)  • 

h(l)+h(2)+h(3)  J  0 

i=l  ' 

1 

•  1 

• 

'  h(l) 

1 

• 

|  h(l)+h(2) 

9 

£  h(i) 

6 

£  h(i) 

1  h(l)+h(2)+h(3) 

| 

1=1  1 

i=l 

1 

Q  =  Q3  =  Diagonal [1,2,4,8,16,32, ..,256]  (A-2) 

The  following  vectors  and  matrices  are  defined: 


*1 


-2 


h(l) 

h(l)+h(2) 

h(l)+h(2)+h(3) 


h(l)+h(2)+..+h(4) 
h(l)+h(2)+..+h(5) 
h(l)+h(2)+. .+h(6) 


-3  -L 


h(l)+h(2)+..+h(7) 

h(l)+h(2)+..+h(8) 

h(l)+h(2)+..+h(9) 


(A- 3) 


(A-4) 


(A-5) 


Qj  =  Diagonal [124]  (A-6) 

Q2  =  Diagonal [  8  16  32  ]  (A-7) 

Q3  =  Diagonal [  64  128  256  ]  (A-8) 
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and  a  piece  of  the  Normal  Equation  can  be  written  as 


/  >  / 
*2  —3 

Q1  0*  0* 

Xj  0  0 

H'QH  = 

0'  x'  x' 

o*  q2  0* 

212  -1  - 

f  r  A 

o*  o*  q3 

—3  -2  21x 

(A-9) 


where 


0  = 


0*  = 


0 
0 
0 

0  0  0 
0  0  0 
0  0  0 


(A- 10) 


(A-l 1) 


Carrying  out  the  matrix  multiplications  indicated  in  Eq  (A-9)  and 
putting  the  result  into  a  summation  form  yields 


H'QH  = 


3 

1  3 

1  3 

EAQA 

k=l 

|k=2^  k-k  i 

!  eXqa-2 

1  k=3 

3 

*  3 

*  3 

kf2^k-lQk-k 

1  EA-lQkA-2 
|  k=3 

3 

1  3 

1  3 

EA-2Qk-k 

k=3 

1  Z  A-2QA-1 
lk“3 

1  .ZA-2QA-2 

|  k=3 

(A-12) 


which  is  symmetric. 

Equivalent  matrices,  containing  all  of  the  information  found  in 
the  original  H  ,  Eq  (22) ,  and  Q  matrices  can  be  formulated  as 


HEQ  -  f2£1»x2*-**2£LJNSMxL  (A- 13) 

QEQ  =  l<Lj  *3.2  •  •  *  ’^L^NSMXL  (A-14) 


where,  for  the  example  problem  Xj,x2  anc*  — 3=— L  are  8iven  by 
Eqs  (A-3)  through  (A-5) ,  and 


*1 


*2 


8 

16 

32 


%  =  %  = 


64 

128 

256 


(A-15) 


(A- 16) 


(A- 17) 


Using  these  equivalent  matrices  (  HEQ  and  QEQ  )  an  inner 
product  relationship  can  be  used  to  produce  [H'QH]  instead  of  the 
straightforward  transpose  and  matrix  multiplications  of  Eqs  (A-l)  and 
(A-2) .  The  first  column  of  the  original  H  matrix  is  formed  using  the 
relationship 

h(i)  -  CFi_1G  (10) 

These  elements  are  put  into  HFC  ,  an  (  NSM*L  )  single  dimension  array. 
HEQ  is  dimensioned  as  an  (  NSMXL  )  matrix  and  equivalenced  to  HFC  . 
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The  FORTRAN  code  used  to  find  iH'QH]  using  HEQ  and  QEQ  is  in 
Appendix  B  within  Subroutine  HBAR. 

At  this  point  the  inverse  of  [H'QH]  is  required.  The  [H'QH] 
matrix  may  be  less  than  full  rank,  requiring  a  generalized  inverse. 

A  library  subroutine,  GMINV,  computes  the  generalized  inverse  through  a 
process  of  column  orchogonalization  using  a  Gram-Schmidt  procedure 
(Ref  13). 

ZEQ  ,  an  (  NSMXL  )  matrix  is  equivalenced  to  Z  ,  an 
(  NSM*L  )  single  dimension  array  containing  the  elements  found  using 
Eq  (21). 


ZEQ  =  lzx,z2>. . ,zL] nsmxl  (A- 18) 

where,  for  the  example  problem 


~yda>  -  ^r 

*1 

yd(2)  -  yzi<2> 

yd(3)  "  yzi(3) 

(A- 19) 

yd(4)  "  yzi(4) 

—2 

yd(5)  -  yzi(5) 

yd(6)  -  yzi(6) 

(A-20) 

”’d<7>  "  yzi(7” 

— 3  e  — L  " 

yd(8)  -  yzi(8) 

yd(9)  "  yzi(S) 

(A-21) 
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An  expression  similar  to  Eq  (A-12)  can  be  built  and  an  inner 
product  relationship  similar  to  the  one  used  to  find  [H'QH]  can  be 
used  to  calculate  [H^Qz]  .  The  FORTRAN  code  used  to  find  [H^Qz] 
using  HEQ  ,  QEQ  and  ZEQ  is  in  Appendix  B  within  Subroutine  CONTROL. 

As  an  indication  of  the  savi  .gs  realized  when  using  the  HEQ  , 
QEQ  and  ZEQ  matrices  versus  H  ,  Q  and  £  ,  consider  the  solution 
of  the  Normal  Equation  for  twenty  step  ahead  prediction  (  L=20  )  and 

ten  smoothing  terms  (  NSM=10  ) .  Matrix  dimensions  are 

^200X10  versus  £HE<^  20X10 

^200X200  versus  [QEQ]  2qx10 

Core  memory  requirements  (in  base  10)  for  these  two  matrices  alone  are 
42000  words  required  for  H  and  Q  versus  400  words  using  HEQ  and 
QEQ  . 
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Appendix  B 


Sample  FORTRAN  Code  Used  in  the  Investigation  of  the  Algorithm 
Implemented  as  a  Regulator 

The  basic  program  in  this  Appendix  follows  the  functional  block 
diagram  of  Figure  1.  Addition  of  a  set  of  logical  switches  allows 
investigation  of  model  mismatch  (Chapter  IV)  and  the  effects  of  noisy 
measurements  and  inputs.  Comment  statements  at  the  front  of  the  program 
explain  the  use  of  these  switches. 

The  only  external  to  the  program  is  a  subroutine  used  to  find 
the  generalized  inverse  of  a  matrix  (GMINV) .  This  subroutine  is  on  a 
functional  library  named  CONTROL,  ID=L  720033,  SN=AFML.  CONTROL  must 
be  attached  and  libraried  before  compilation  of  the  program. 

The  output  of  this  particular  version  of  the  algorithm  is 
Figures  28  through  31  in  Chapter  III. 
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Appendix  C 


Graphical  Results  of  Model  Mismatch  Simulation 


Appendix  C  contains  the  graphical  results  corresponding  to  the 
cases  listed  in  Table  II  and  Table  III  in  Chapter  IV. 
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SAMPLED  OUTPUT 
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oamtieoci  re  m 
>  hmm  reue  re  ni 

MR>4  L*7  TOO.  in 4 

•*iTc«c4«T-r.T.r.r 


Fig  C-l.  System  Output  for  Truth  Model  (A)  and  10%  Delta  1  Model 
to  Zero- Input  Response  Calculation,  TC=0.075 
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Fig  C-2.  Controls  Applied  for  Truth  Model  to  Zero-Input 
Response  Calculation,  TC=0.075 


SAMPLED  OUTPUT 


CONTROL 


o 

o 


Fig 


C-7.  System  Output  for  Truth  Model  (A)  and  10%  Delta  1  Model  to 
Zero-Input  Response  Calculation,  TC=0.124 


Fig  C-8 


Controls  Applied  for  Truth  Model  to  Zero-Input 
Response  Calculation,  TC=0.124 


Fig  C-9.  Controls  Applied  for  10%  Delta  1  Model  to  Zero-Input 
Response  Calculation,  TC-0.124 
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SAMPLED  OUTPUT 


Fig  C-10.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model 
to  Zero-Input  Response  Calculation,  TC=0.075 


Fig  C-ll.  Controls  Applied  for  Truth  Model  to  Zero-Input 
Response  Calculation,  TC=0.075 


Fig  C-13.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model 
to  Zero-Input  Response  Calculation,  TC=0.088 
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TIME  IN  SECONDS 


Fig  C-16.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model 
to  Zero-Input  Response  Calculation,  TC=0.124 


Fig  C-18.  Controls  Applied  for  20%  Delta  2  Model  to  Zero-Input 
Response  Calculation,  TC=0.124 
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SAMPLED  OUTPUT 

H.QO  0.00  4.00  8.00  12.00  16.00  20.00  24.00 


Fig  C-19.  System  Output  for  Truth  Model  (A)  and  30%  Delta  3  Model  to 
Zero-Input  Response  Calculation,  TC=0.075 
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Fig  C-20.  Controls  Applied  for  Truth  Model  to  Zero-Input 
Response  Calculation,  TC=0.075 


Fig  C-21.  Controls  Applied  for  30%  Delta  3  Model  to  Zero-Input 
Response  Calculation,  TC=0.075 
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SPilPLED  OUTPUT 

0.00  4.00  8.00  12.00  16.00  20.00  24.00 
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SHITCHEfliT-F.T.F.F 


Fig  C-22.  System  Output  for  Truth  Model  (A)  and  30%  Delta  3  Model  to 
Zero-Input  Calculation,  TC=0.088 
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Fig  C-23.  Controls  Applied  for  Truth  Model  to  Zero-Input 
Response  Calculation,  TC=0.088 


Fig  C-24.  Controls  Applied  for  30%  Delta  3  Model  to  Zero-Input 
Response  Calculation,  TC=0.088 


125 


SAMPLED  OUTPUT 

*-4  •  00  0.00  4.00  8.00  12.00  16.00  20.00  24.00 


.00  0.40  0.80  1.20  1.60 


TIME  IN  SECONDS 


Fig  C-25.  System  Output  for  Truth  Model  (A)  and  10%  Delta  1  Model 
to  H  Calculation,  TC=0.075 
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2.  BO 

TIME  IN  SECONDS 

Fig  C-26.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.075 
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0.40 

0.  BO  1.20 

1.60  2.00 

2.40 

2.80 
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SECONDS 

Fig  C-27 .  Controls  Applied  for  10%  Delta  1  Model  to  H 
Calculation,  TC=0.075 
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SAMPLED  OUTPUT 

r4.00  0.00  4.00  8.00  12.00  16.00  20.00  24.00 


Fig  C-28.  System  Output  for  Truth  Model  (A)  and  10%  Delta  1  Model  to 
H  Calculation,  TC=0.088 
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Fig  C-29.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.088 


Fig  C-30.  Controls  Applied  for  10%  Delta  1  Model  to  H 
Calculation,  TC  =  0.088 
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SAMPLED  OUTPUT 

r4.00  0.00  4.00  8.00  12.00  16.00  20.00  24.00 


.00  0.80  1.60  2.40  3.20 

TIME  IN  SECONDS 


Fig  C-31.  System  Output  for  Truth  Model  (A)  and  10%  Delta  1  Model  to 
H  Calculation,  TC=0.124 
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Fig  C-32.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC  =  0.124 


Fig  C-33.  Controls  Applied  for  10%  Delta  1  Model  to  H 
Calculation,  TC=0.124 
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SAMPLED  OUTPUT 

r6.00  -2.00  2.00  6.00  10.00  14.00  18.00  22.00 


Fig  C-34.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model  to 
H  Calculation,  TC=0.075 
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Fig  C-35.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.075 


Fig  C-36.  Controls  Applied  for  20%  Delta  2  Model  to  H 
Calculation,  TC=0.075 
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Fig  C-37.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model  to 
H  Calculation,  TC=0.088 
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QQ  0.40  0.80  1.20  1.60  2.00  2.40  2.80 

TIME  IN  SECONDS 


Fig  C-38.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.088 
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2.80 

TIME  IN  SECONDS 

Fig  C-39.  Controls  Applied  for  20%  Delta  2  Model  to  H 
Calculation,  TC=0.088 


SAMPLED  OUTPUT 
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.00  0.80  1.60  2.40  3.20 
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Fig  C-40.  System  Output  for  Truth  Model  (A)  and  20%  Delta  2  Model 
to  H  Calculation,  TC=0.124 
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I  Fig  C-43.  System  Output  for  Truth  Model  (A)  and  30%  Delta  3  Model  to 

|  H  Calculation,  TC=0.075 
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Fig  C-44.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.075 
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Fig  C-45.  Controls  Applied  for  30%  Delta  3  Model  to  H 
Calculation,  TC-0.075 
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SAMPLED  OUTPUT 

.-7.00  -3.00  3.00  5.00  9.00  33.00  37.00  21.00 


Fig  C-46.  System  Output  for  Truth  Model  (A)  and  30%  Delta  3  Model  to 
H  Calculation,  TC=0.088 
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*  HRRK8  TRUE  TO  HQflR 
H5h=4  L«7  TC* 0.088 
MITCHCS*  T.  T-F.F.F 


Fig  C-47.  Controls  Applied  for  Truth  Model  to  H 
Calculation ,  TC=0 . 088 


SAMPLED  OUTPUT 

.-10.00  -5.00  0.00  5.00  10.00  15.00  20.00  25.00 


Fig 


1-49.  System  Output  for  Truth  Model  (A)  and  30%  Delta  3  Model  to 
H  Calculation,  TC-0.124 
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Fig  C-50.  Controls  Applied  for  Truth  Model  to  H 
Calculation,  TC=0.124 


Fig  C-51.  Controls  Applied  for  30%  Delta  3  Model  to  H 
Calculation,  TC=0.124 
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Appendix  D 


FORTRAN  Code  for  Terrain  Following  Test  Cases  of  Chapter  V 

This  Appendix  contains  the  code  used  in  the  terrain  following 
Test  Case  1  of  Chapter  V.  Two  libraries  must  be  attached  prior  to 
compilation: 

1.  1MSL,  ID  =  Library,  SN  =  ASD  for  a  sorting  routine. 

2.  CONTROL,  ID  =  L720033,  SN  =  AFML  for  a  matrix  inversion 

routine. 

Figures  50  through  53  are  outputs  from  this  program. 
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